% The code is based on the following paper:
%   Chen, T. Y., Lin, Y. L., and Tzeng, L. Y., forthcoming.
%   Estimating probability weighting functions through option pricing bounds. Review of Asset Pricing Studies.
%
% Copyright: Tzu-Ying Chen, Yo-Lan Lin, Larry Y. Tzeng
% Date: January 30, 2024

clear

%% Setting           
BV_Target = {'UB'};                                                        

%% Optimal Risk Preferences under Rank-Dependent Expected Utility
FileName = ['Summary_OPT_AllCR_AllPeriod_OP_UB'];
load(FileName, ...
     'Date_Monthly', ...
     'Index_CR', 'Index_OPT', 'OPT_TAB');
clear FileName 

OPT = nan(size(Date_Monthly, 1), length(Index_OPT));                       
for d = 1:size(Date_Monthly, 1)
    OPT(d, :) = OPT_TAB(Index_CR, Index_OPT, d);                                 
end 
clear Index_CR Index_OPT OPT_TAB
clear d

%% Four Specific Cases 
AllPeriod = [[19990317 19990416]; ...
	         [20030219 20030321]; ...
	         [20111019 20111118]; ...
	         [20090121 20090220]];      
[~, Index_Location] = ismember(AllPeriod, Date_Monthly(:, [1 2]), 'rows');
Index_Location(Index_Location==0) = [];                                    
clear Date_Monthly

OPT = OPT(Index_Location, :);                                              
clear Index_Location

AllPeriod_Text = cell(size(AllPeriod, 1), 1);                              
for d = 1:size(AllPeriod, 1)
    if d < size(AllPeriod, 1)
        AllPeriod_Text(d, 1) = cellstr(char(['\qquad ' datestr(datenum(num2str(AllPeriod(d, 1)), 'yyyymmdd'), 'mmm dd, yyyy') ' (' num2str(25 * d) 'th percentile of $\alpha_{t}^{\ast}$)']));
    else
        AllPeriod_Text(d, 1) = cellstr(char(['\qquad \quad ' datestr(datenum(num2str(AllPeriod(d, 1)), 'yyyymmdd'), 'mmm dd, yyyy') ' (maximum of $\alpha_{t}^{\ast}$)']));
    end
end
clear d

%% Plot Figure: Probability Weighting Function (Slope)
z = 0:0.001:1;
eps = 10^(- 6);

% Plot Figure
h_PWF = figure;
for d = 1:size(AllPeriod, 1)
    FileName = ['Summary_OP_Grid_' num2str(AllPeriod(d, 1))];
    load(FileName, ...
         'Record');
    clear FileName

    Index = find((Record(:, 1)==1) & ...
                 (Record(:, 2)==1));
    VIORate = 1 - Record(Index, end - 1);                                                             
    clear Index Record
    
    % Around Low Return (Bisection Method)
    x_L = min(z) + eps;
    x_U = median(z) - eps;
    x_M = 0.5 * (x_L + x_U);   
    for n = 1:10000
        x = [x_L x_M x_U];
        y = P1998_D1(x, OPT(d, 1), OPT(d, 2)) - P1998_D1(x, 1, 1);
                    
        if (sign(y(1)) * sign(y(2))) < 0
            x_U = x_M;                                                     
        elseif (sign(y(2)) * sign(y(3))) < 0
            x_L = x_M;                                                     
        else
            x_sol = nan;
            break            
        end
        x_M = 0.5 * (x_L + x_U);                                           
            
        Error = P1998_D1(x_M, OPT(d, 1), OPT(d, 2)) - P1998_D1(x_M, 1, 1);
        if abs(Error) < eps
            x_sol = x_M;
            break
        else
            x_sol = nan;
        end
    end
    Point_LowRET = x_sol;
    clear x_L x_U x_M
    clear x y Error   
    clear x_sol

    % Around High Return (Bisection Method)
    x_L = median(z) + eps;
    x_U = max(z) - eps;
    x_M = 0.5 * (x_L + x_U);   
    for n = 1:10000
        x = [x_L x_M x_U];
        y = P1998_D1(x, OPT(d, 1), OPT(d, 2)) - P1998_D1(x, 1, 1);
                    
        if (sign(y(1)) * sign(y(2))) < 0
            x_U = x_M;                                                     
        elseif (sign(y(2)) * sign(y(3))) < 0
            x_L = x_M;                                                     
        else
            x_sol = nan;
            break            
        end
        x_M = 0.5 * (x_L + x_U);                                           

        Error = P1998_D1(x_M, OPT(d, 1), OPT(d, 2)) - P1998_D1(x_M, 1, 1);
        if abs(Error) < eps
            x_sol = x_M;
            break
        else
            x_sol = nan;
        end
    end
    Point_HighRET = x_sol;
    clear x_L x_U x_M
    clear x y Error   
    clear x_sol
    
    % Plot Figure
    subplot(2, 2, d)
    plot(z, P1998_D1(z, OPT(d, 1), OPT(d, 2)), ...
         'LineStyle', '-', ...
         'LineWidth', 2, ...
         'Color', 'r')
    hold on

    plot(z, P1998_D1(z, 1, 1), ...
         'LineStyle', ':', ...
         'LineWidth', 2, ...
         'Color', 'k')
    hold on
    grid on

    text(Point_LowRET, 0.25 + P1998_D1(Point_LowRET, 1, 1), ...
         ['$\uparrow \underline{p}$ = ' num2str(roundn(Point_LowRET, - 4))], ...
         'HorizontalAlignment', 'Left', ...
         'VerticalAlignment', 'BaseLine', ...
         'FontSize', 15, ...
         'FontName', 'Times New Roman', ...
         'FontWeight', 'Bold', ...
         'LineWidth', 2, ...
         'Color', 'k', ...
         'Interpreter', 'Latex')          
    hold on    
    
    text(Point_HighRET, 0.25 + P1998_D1(Point_HighRET, 1, 1), ...
         ['$\uparrow \overline{p}$ = ' num2str(roundn(Point_HighRET, - 4))], ...
         'HorizontalAlignment', 'Left', ...
         'VerticalAlignment', 'BaseLine', ...
         'FontSize', 15, ...
         'FontName', 'Times New Roman', ...
         'FontWeight', 'Bold', ...
         'LineWidth', 2, ...
         'Color', 'k', ...
         'Interpreter', 'Latex')          
    hold on  
    
    h = legend(['($\alpha^\ast$ , $\beta^\ast$) = (' num2str(OPT(d, 1)) ' , ' num2str(OPT(d, 2)) ')']);
    set(h, 'FontSize', 15, ...
           'FontName', 'Times New Roman', ...
           'FontWeight', 'Bold', ...
           'Location', 'NorthWest', ...
           'Interpreter', 'Latex', ...
           'Box', 'Off')
    clear h
    
    h1 = xlabel('$p$');  
    h2 = ylabel('slope of $v(p)$');   
    h3 = title([AllPeriod_Text(d, :); '(percentage of violations under EU = ' num2str(roundn(100 * VIORate, - 2)) '\%)']);     
    set([h1 h2 h3], 'FontSize', 15, ...
                    'FontName', 'Times New Roman', ...
                    'FontWeight', 'Bold', ...
                    'Interpreter', 'Latex')
    clear h1 h2 h3    
  
    set(gca, 'XLim', [0.005 0.995], ...
             'XTick', [0.005 0.25 0.5 0.75 0.995], ...
             'FontSize', 15, ...
             'FontName', 'Times New Roman')
    set(gca, 'YLim', [0 7], ...
             'YTick', 0:1:7, ...
             'FontSize', 15, ...
             'FontName', 'Times New Roman')    
    
    set(gca, 'Layer', 'Top')
    clear VIORate Point_LowRET Point_HighRET
end
clear d n

X = 29.7;                                                                  
Y = 21.0;                                                                  
XMargin = 1;                                                               
YMargin = 1;                                                               
XSize = X - 2 * XMargin;                                                   
YSize = Y - 2 * YMargin;                                                   
set(h_PWF, 'PaperUnits', 'Centimeters')
set(h_PWF, 'PaperSize', [X Y])
set(h_PWF, 'PaperPosition', [XMargin YMargin XSize YSize])
set(h_PWF, 'PaperOrientation', 'Portrait')
clear AllPeriod AllPeriod_Text AllCR_Target h_PWF X Y XMargin YMargin XSize YSize
clear z eps

%% Output
FileName = ['Figure_PWF_Slope_UB'];
saveas(gcf, [FileName '.eps'], 'epsc')
print(FileName, '-dpdf', '-fillpage')
winopen([FileName '.pdf'])
clear FileName